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We consider the non-equilibrium dynamics of the interacting Lieb-Liniger gas after instantaneously 
switching the interactions off. The subsequent time evolution of the space- and time-dependent 
correlation functions is computed exactly. Different relaxation behavior is observed for different 
correlation functions. The long time average is compared with the predictions of several statistical 
ensembles. The generalized Gibbs ensemble restricted to a fixed number of particles is shown to 
give correct results at large times for all length scales. 



With the recent advances in ultracold atomic gases 
it is now possible to realize isolated quantum systems 
with long coherence times. These experiments are ideal 
to study non-equilibrium physics [IH3]. The outcome of 
these experiments initiated the ongoing debate concern- 
ing the circumstances under which an isolated quantum 
system will thermalize, which has led to many theoreti- 
cal studies on non-equilibrium dynamics of quantum sys- 
tems. Most of these studies considered a so-called quan- 
tum quench: an isolated system is prepared in the ground 
state of a Hamiltonian H(c), where c is some control- 
lable parameter, like an interaction strength. At to = 
the parameter c is instantaneously switched to a differ- 
ent value resulting in unitary time evolution governed by 
the new Hamiltonian. Considerable progress has been 
made (for a recent review see [4] ), however there are 
still many open questions. Under what conditions can a 
system equilibrate and what is the relevant mechanism? 
What is the role played by integrability or its absence? 
Do different kinds of correlation functions show different 
relaxation behavior? And what is the importance of the 
initial state? 

The difficulty with answering these questions is that 
one must rely on case by case studies, thus, one is faced 
with the problem of distinguishing universal behavior 
from case-specific results. Moreover, most methods used 
so far suffer from the fact that they are either valid for 
short times or for large ones. What we propose here is to 
study a specific example where we can study the behavior 
of correlation functions at all time scales, using an exact 
method. We can, therefore, not only determine what the 
equilibrium state is, if there is one, but also how it is 
reached. Furthermore, the method we use allows us to 
consider various types of correlation functions, thereby 
obtaining a more complete picture for this case. 

The model we consider is the Lieb-Liniger model, 
which has experimentally been realized in various circum- 
stances [2JG2HJI6]. We bring the model out of equilibrium 
by instantaneously turning off the interactions; this is a 
special case of an interaction quench. For short times one 
can think of this as a simulation of experiments where an 
ultracold Bose gas is released from a tight transverse con- 
finement p3[B]. 



Several studies on the non-equilibrium dynamics of the 
Lieb-Liniger model have appeared before. The effect of 
instantaneously turning on the interactions was inves- 
tigated in 9 using the Bethe Ansatz in combination 
with a Monte-Carlo sampling technique, and in |10j using 
the numerical time-evolving block decimation algorithm. 
The expansion starting from a regular array was studied 
in [llj . where the properties of this specific initial state 
were exploited using the coordinate Bethe Ansatz. In[T2] 
the possibility of studying quenches in the Lieb-Liniger 
model using an integrable field theory is discussed. 

The outline of this article is as follows. In section U 
we describe the setup and explain the methods to com- 
pute correlation functions after the quench. This is fol- 
lowed by section |TT] where we give the results for the time 
evolution of the non-local pair correlation and the auto- 
correlation. We also introduce a new type of correlation 
function that is a measure for the correlation between the 
pre- and post-quench states, which we will call quench- 
straddling correlations. In section |III| we compare the 
long time average of the correlation functions with the 
predictions of various statistical ensembles. First we dis- 
cuss why the Generalized Gibbs (GGE) ensemble fails in 
this case. Considering the generalized canonical ensem- 
ble (GCE), that is the GGE but with the total number of 
particles explicitly fixed, we show that it yields the cor- 
rect results at all length scales. We end with a discussion 
and conclusions. 



I. SETUP AND METHOD 

The model we discuss in this paper is the one- 
dimensional Bose gas which is described by the Lieb- 
Liniger Hamiltonian |13) 

H(c)= [ dx{d x ¥(x)d x V{x) + c^(x)^{x)^(x)^(x)} 
Jo 

(1) 

with L the length of the system and c the interaction 
strength. For simplicity we will work with periodic 
boundary conditions. We will only consider the repul- 
sive case c > 0. Two limiting regimes of the Lieb-Liniger 
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model are: noninteracting Bosons (c = 0), and the so- 
called Tonks-Girardeau regime (c = oo) [H] where the 
bosons effectively behave as free fermions. In order to 
study the non-equilibrium dynamics, we quantum quench 
the system [T3 [H]. We do this by first preparing the 
system in the ground state for c > and then instan- 
taneously turning off the interactions: H(c) H(0) at 
t = 0. 

The equal-time 2-point function (* t (a;)*(0)) will not 
evolve in time for this specific quench. The Fourier- 
transform of this 2-point function consists of computing 
(W k ^ q ). Since the states before and after the quench are 
translationally invariant, one only needs to compute the 
diagonal components (q = k). However, when c = 0, the 
momentum occupation operators ^^j. are constants of 
motion, hence the 2-point function will not be affected by 
this quench. For higher-point correlations or dynamical 
ones this will not be the case. As an illustration we will 
consider the measurable non-local pair correlation 



92 (x) = 



(* t (a;)* t (0)^(a;)*(0)) 
(*t(Q)^(0)) 2 



(2) 



The non-local pair correlation for finite size systems 
has been computed in [17 using exact methods for the 
ground state. The case of finite temperature was stud- 
ied in [TS] using perturbative techniques. Without loss 
of generality, we will fix in the following the density at 
unity (*t(o)*(0)) = N/L = 1. 

In order to study the time evolution after the interac- 
tions are turned off, it is useful to write the final Hamil- 
tonian as H(0) = Wfe^t^t, with dispersion relation 
ujk = k 2 . The time evolution of the field operators read- 
ily follows: e ^(o)t^t e -^(o)t = $t e *fc 2 t. So instead of 
acting with the time evolution operator on the initial 
state, we can act with it on the operators. The upshot 
of this is that we can avoid decomposing the initial state 
in terms of eigenstates of the final Hamiltonian. This 
decomposition typically involves a number of states that 
grows exponentially with the total number of particles. 

Using the Fourier-transform of the field operators 
*t( a 



^ fc e %kx W k the Heisenberg picture allows us 
to express the time evolution of the dynamical correlation 
function after the quench as 



g 2 (x,t 1: t 2 ) 
1 



E 

fei,fc 2 ,fc3 



= (4>\&(x, t 2 )* t (O,i 1 )*(a;,t 2 )*(O,t 1 )|0) 



/(x,ta s *i;{*i}) = MaJ-2A?a3ti + (fri+fcs)(ta-ti)) ( 3 ) 

where kij = ki — kj and we used k\ + k 2 = k 3 + k± be- 
cause of translational invariance. The initial state 
is the ground state of the fully interacting Licb-Linigcr 
model (c > 0). The quench takes place at t = and we 
evaluate the correlation function at times < t\ < t 2 . 



In equilibrium (at t 2 = t% — to), this is the normal or- 
dered 2-point function of current ^(x)^(x) operators 
and can efficiently be obtained using the techniques de- 
scribed in [17] . In order to study the time evolution, we 
need to compute a four-point function of field operators 
^t. This four-point function of the initial state \<f>) can 
be computed in the framework of the Algebraic Bethe 
Ansatz [19] by generalizing the methods of [20] for their 
computation of two-point functions. We will briefly dis- 
cuss how the computations are done without going in too 
much detail. First, we express the correlation function in 
terms of matrix elements. By inserting a resolution of 
the identity operator between every pair of adjacent field 
operators, one can write the four-point function in ^ as 



J2 (0|* t JrM)(m|<K> 

1,"2,"3 

(n 2 \y k3 \n 3 )(n 3 \y kl ^ 2 _ k3 \cj)) (4) 



in terms of the matrix elements (n|\I/jj.|m) of the field op- 
erator. Here \(f>) is the initial state for N particles. The 
states \n 2 ) and \n 3 } are intermediate states with 

N — l, N—2 and N — l particles respectively. In order 
to compute the matrix element (nl^ljm) one first solves 
the Bethe equations for the states |n) and |m) which 
results in a set of so-called rapidities for both states. 
The matrix elements can then be evaluated by comput- 
ing the determinant of a matrix whose entries are ra- 
tional functions of the rapidities of the two eigenstates 
involved. The explicit expression for the matrix element 
can be found in [5T] as a sum of determinants and in [5U] 
(see eq. 41) in terms of a single determinant. What re- 
mains to be performed are the actual summations over 
Til, n 2 and n 3 . Starting from the initial state \<f>) the 
Fock space of intermediate states \ni) is scanned by nav- 
igating through choices of sets of quantum numbers. For 
each individual intermediate state, the Bethe equations 
are solved, and the matrix element is computed. The 
search for important states is done via the ABACUS 
method [55], and is close to optimal. In order to verify 
the accuracy of our results, we keep track of the sum-rule 

l(0l*fcl n i)| 2 = N/L. The sum over m is truncated 
after a desired precision is achieved. Next, for every state 
|ni), we repeat the process by performing a summation 
over n 2 , for which we can also compute a sum-rule. We 
do not need to perform the summation over n 3 explic- 
itly because we can use the hermitian conjugate of the 
matrix element (0|^ / jj ;i ^ , fc 2 |"-2), thereby constructing the 
full four-point function. For a system of N — 20 particles 
and large interaction c = 1000, which is the most difficult 
case, a number of intermediate states \n 2 ) of the order of 
10 s are needed in order to saturate the sum-rule at 99.6%. 
Once the four-point function is obtained, dynamical cor- 
relation functions ^ can be computed straightforwardly. 
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II. RESULTS 

In this section, we study the dynamical correlation 
function .92(^,^1,^2) ^ after the quench, for all times. 
We present the results for a system of size L = N = 20 
starting from the ground state for various c > 0. For 
clarity, we specialize to the cases of the non-local pair 
correlation g2(x,t) — g2(x,t,t) and the auto-correlation 
02 £2) = 52(0, ti,t 2 ). We conclude this section by con- 
sidering what we will call the quench-straddling correla- 
tions. 



A. Non-local pair correlation 

1. Short Times 

First, we consider the non-local pair correlation. We 
plot gi (x, t) for various times t after the quench as a func- 
tion of x in fig. [I] For clarity we also plot how g 2 (x,t) 
for various fixed values of x, develops in time in fig. [2] 
We display here only the results for c = 1000 because 
for smaller c the behavior is roughly the same but less 
pronounced. In equilibrium (t — 0) we see that g 2 (x) van- 
ishes for small result of destructive interference 
due to the fermionic character of the Tonks- Girardeau 
gas (this is also called anti-bunching) . For large x, cor- 
relations decay, no destructive interference takes place 
and g2{x) — > 1. If we now focus on the behavior for 
x = we see that because of loss of coherence the cor- 
relation builds up. Eventually, for large times there will 
be constructive interference ( g 2 (x) > 1), which is called 
bunching, and is typical for free bosons (see fig. [2]). For 
small x and small t we also expect that correlations will 
grow because of loss of coherence, however we have to 
take into account that g 2 (x, t)dx = 1 — 1/N is a con- 
served quantity which results in non-monotonic behavior 
for x > at small times as can be clearly seen in fig. [2j 
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FIG. 1. g2(x,t) as a function of t for various x. N,L = 20 
and c = 1000. 




FIG. 2. g2{x,t) as function of x (solid lines) for various 
times together with asymptotic predictions (dashed lines), for 
N, L = 20 and c = 1000. 
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FIG. 3. \g2{x,t) — g2(x, 0)|, together with x — v g t (dashed 
line) for c = 1000 and N,L = 20. 



In order to study how the correlations propagate 
through the system we plot \g2(x, t) — g2(x, 0)| in fig. [3] 
There is no light-cone effect because of the absence of a 
maximal velocity. As reference we have plotted x = v g t 
(dashed line) where the group velocity is computed as 



2E* 



take place after t • 



) . We see that the dominant changes 
x/v g . 



2. Finite size effects 

From a finite size study we can conclude that the be- 
havior for short times, discussed in the previous sec- 
tion, is independent of the system size. If we now 
turn to large times finite size effects become more and 
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FIG. 4. Complete and partial revival of g2(x = 0, t) compared 
with the fits for a system of finite size and infinite size, while 
keeping the density fixed (N/L — 1). 



more pronounced. Because of the simple dispersion 
relation u(k n ) = (2irn/L) 2 exact revival occurs at 
i rov = L 2 /(Att). Partial revivals of decreasing strengths 
also occur at higher harmonics of the revival time: 
i rov /2, i rcv /3, i rov /4, . . .. Due to the revival, 32(0, t) gets 
suppressed for large times, as can be seen in fig. [4j In 
the thermodynamic limit there will be no suppression 
and one can easily show that 52(0) = 2 for free bosons in 
equilibrium. 



3. Asymptotic behavior 

In order to understand the intermediate time behavior, 
which is after the non-monotonic behavior and before 
finite size effects start to dominate, we fit the results for 
small x with the following test function 



g 2 (x,t;c) = a(x\c) 



b(x; c) 



|i rev sin(7rf/t rev )|«( c ) 



(5) 



We motivate this function as follows. The first part 
a(x] c) is the stationary part. For the time-dependent 
part we expect algebraic decay, since the initial state is 
quantum critical. The exact revival, due to the finite 
size, is taken into account by means of the sine-function. 
The fits for small x and t are compared with the exact 
result in fig. [2j The data suggests that a(x, c) is a de- 
creasing function of x which agrees with the result for 
the long time average which we will discuss in section 
III We find that a(c) is a decreasing function of c and 
a ~ 1/2 for c = 1000. However, since we have only a 
limited range where we can fit the data and a is small, 
no firm predictions can be made. After comparing fits 
for various system sizes, our data seem to be consistent 
with the dependence on L entering only via 
agrees with the intuition that for small x and small times 
t after the quench finite size effects are irrelevant. Assum- 
ing this is indeed the case, one can from the finite size 



data make predictions for the thermodynamic limit by 
sending t rev — » 00 in ^ while keeping the other param- 
eters fixed. The stationary part a(x; c) would then be 
the large time limit in the thermodynamic limit, since 
the time-dependent part now vanishes for large t, which 
is not the case for finite size systems. For free bosons in 
equilibrium one can easily show that 32(0) = 2, which 
is compatible with what we find from the fitting data: 
a(0,c) ~ 2. In fig. Hwe compare the exact results for 
52 (0, t) with the fit both for finite size and our pre- 
dictions for the thermodynamic limit. 

Since the post-quench state is completely described by 
the correlations on the initial ground state, one can try 
to predict asymptotic behavior, such as the exponent in 
([5]), using low-energy effective theories like bosonization. 
If we write ^ as 



92(x,t) = ±J2 f ^ dzetkZ 



($ , (i-2fct)$ t (z)*(z+i-2H)$(0)), (6) 

one can see that because of the integral over z one 
needs the correlation function at all length scales. Hence 
whether low-energy descriptions can be used or not re- 
mains an open question. 



B. Auto-correlation 

The auto-correlation 32(^1,^2) is plotted in fig. [5] as 
a function of t 2 after various times t\ after the quench. 
As for the non-local pair correlation, there is an anti- 
bunching bunching transition. However, in contrast to 
the evolution of the non-local pair correlation, the auto- 
correlation does increase monotonically in time after the 
quench, as can be seen in fig. [5j This is consistent with 
the fact that the integrated auto-correlation is not con- 
served in time, in contrast with the non-local pair corre- 
lations. 



C. Quench-straddling correlations 

The versatility of our method allows the computation 
of various different dynamical correlation functions. As 
an example we will study a correlation function which is 
a measure of how the pre- and post-quench states are cor- 
related, what we will call a quench-straddling correlation 
function. We consider the situation where the quench 
takes place at time to = 0. We then compute a two-point 
function where we evaluate 'F(O) before the quench at 
time t- < to and ^(x) after the quench, at time i+ > to 

3st r addic(x,t_,i + ) = (d ) \e m ^¥(x)e- iH ^ 

e -ii?( C )t_^^ e iH(c)t_|^_ (7) 



To simplify this expression we first use the Fourier- 
transform of the field operators ^(x) — -/fJ2k e 
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FIG. 5. The auto-correlation 52(0, £2; ii) for various times 
ti = 0.00,0.05,0.10,0.14,0.24,0.53 as a function of t 2 for a 
system with N, L = 20 and c = 1000. 
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FIG. 6. The quench-straddling correlation function 
Re{ 5 

straddle (0,t-,t+)} as a function of t+ — t-. The vari- 
ous graphs (from top to bottom) correspond to the quench 
times to — t- — 0.0, 0.2, 0.4, 0.7, 1.4, 00, indicated by the ver- 
tical lines. The system parameters are N, L — 80 and c = 5. 



As before, we can then handle the time-evolution of 
the post-quench Hamiltonian H(Q) using the Heisenberg- 
picture. For the time evolution of the pre-quench 
Hamiltonian H{c) we use the Schrodinger picture; this 
is achieved by inserting a resolution of the identity 
S n l n )( n li m t erms °f eigenstates |n) of the pre-quench 
Hamiltonian between e~ lH ^ t+ and e^ lH ^ 1 - . By 
pulling out all the different phases and using translation 
invariance we arrive at 



Astraddle (X, t-,t + ) 

1 



E i(E n -E )t-+ik^t + +ik n x 



*t \n) 



(8) 
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FIG. 7. The long time average g2(x) for various c for a system 
with N, L — 20. 



Here Eq is the energy of the initial state \<j>) and E n , k n 
are the energy and momentum of the intermediate state 
\n) respectively. In fig. [6] we plot Re{g str addic(0, t+)} 
for various t_ as a function of t = t + . A first observation 
is that in the case of no quench, t+ — to, the correla- 
tion decays much faster than the case after the quench 
t- = to- This can be understood from that fact that, 
for a sufficiently large interaction strength c, the disper- 
sion E(k) is linear for small k. Therefore, the majority 
of the states \n) in |7]) have an energy E(k) > k 2 . When 
considering the general quench times to — t_, a striking 
observation we can make is that just after the quench the 
correlation function suddenly seems completely relaxed. 
The explanation is that for a quench-straddling correla- 
tion dephasing takes place both via the energy and the 
momenta, which can be considered as orthogonal direc- 
tions. This speeds up the relaxation significantly for a 
short period right after the quench. The only coherence 
that is left is for very small k, resulting in an extremely 
slow relaxation compared to the case if there were no 
quench. 



III. STATISTICAL ENSEMBLES 

In this section we discuss whether the large time be- 
havior of correlation functions after the quench can be 
described by a statistical ensemble. For simplicity we re- 
strict ourselves in this section to the analyses of g2(x,t). 
Due to the finiteness of the system, no actual relaxation 
occurs. However, for most of the time, g2(x,t) oscillates 
around the same mean value, as is seen in fig. |4j There- 
fore, it is still useful to consider the long time average 
(LTA) of correlation functions 



g 2 (x;c)= lim -/ g 2 (x,t;c)dt. (9) 

J ->UTJ 1 Jo 

Note that because of the exact revival we can take 
Tlta = trev ■ The result of the time average is presented 
in fig. [7J It is well-known that quantum mtegrable mod- 
els do not always relax to a state of thermal equilibrium. 
As an alternative the Generalized Gibbs ensemble (GGE) 



G 



[2"3"] has been proposed. The GGE is based on the idea of 
Jaynes [24] to construct a statistical ensemble with max- 
imal entropy subject to constraints for the expectation 
values of the conserved charges. Expectation values of 
observables in the GGE are then computed via 



(0(x)) 



GGE 



= Tr 



{0(a;)e-£"' 9 " Q "} /Z, 



GGE 



(10) 



where Zqge — Tr {e~ ^« ^"Q™ } is the generalized par- 
tition function. The Langrange multipliers fa corre- 
sponding to the conserved charges Q n are fixed via 
the initial conditions (4>\Q n \4>) = (Qh)gge- So far, 
the GGE has only been used for models which are ef- 
fectively free, namely the Hamiltonian can be written 
as H = J\ uik^f\^f k , with the momentum occupation 

numbers ^\.^ k as the obvious choice for the conserved 
charges. Many successful examples exist [TBI |2"3"I |2"5H3"T] . 
Usually, the GGE is implemented as a grand canonical 
ensemble, which simplifies the results because the parti- 
tion function can be written as a product of single states. 
For example, the partition function can be written as 



J GGE 



IK 1 *- 



-A\±i 



(11) 



where the + and — signs corresponds to fermions and 
bosons respectively. The Lagrange multipliers are deter- 
= l/(e /3fc ± 1), (+ fermions and 



mined via 



— bosons). However, in this procedure the correlation 
between different conserved charges are not kept; that is, 

<***fc*J*g>GGB = (*!* fc >GGB<*J* g )GGE for k + q. 

Despite these shortcomings, the GGE has been successful 
in many cases, although exceptions are known, for exam- 
ple when translation invariance is broken [33] [33J . Let 
us now turn to the case of the construction of the GGE 
for the LTA of g 2 (x). First we explicitly write down the 
LTA of <?2 ( x ) using (|} and {9} 



ki^k 2 

■ N(N-l) 
^ L 2 



(12) 



One can see that the LTA for the g 2 (x) is explicitly writ- 
ten as a sum over expectation values of products of the 
conserved charges. From this we can see that the GGE 
is not applicable in this case, since the initial state is the 
ground state of an interacting system and Wick's theo- 
rem does not apply here. Instead we considered the gen- 
eralized canonical ensemble (GCE) by keeping the total 
number of particles fixed. 



A. The generalized canonical ensemble 



Consider a system where the Hamiltonian can be diag- 



We 



onalized in terms of free particles H = ^ fc UkWJ& k . 
impose no restrictions on the dispersion relation ujk and 



the operators can have either fermionic or bosonic 
commutation relations. The partition function for such 
a system in the canonical ensemble (CE) can be written 

as 



(13) 



no— ni 



=0 k 



where n max = 1,N for fermions and bosons respectively. 
The presence of Sj2-nj,N makes it very complicated to 
evaluate the sums directly. Fortunately, one can compute 
the partition function Zn for a system of N particles via 
a well known recursion relation (see for instance |34H36j ) 



1 



N 



71=1 



Zn- 



(14) 



with Z n — for n < 0, Zq = 1 and z n = J2k ex P[ — n P UJ k\- 
The minus and plus signs correspond to fermions and 
bosons respectively. For a given ft and N the parti- 
tion function can be evaluated numerically. The com- 
putational complexity is of the order of n cu toff^V 2 , where 
^cutoff is the number of one-particle states considered in 
the computation of z n . One can easily generalize the 
computation of the canonical ensemble to what we will 
call the generalized canonical ensemble (GCE) by making 
the replacement f3uj)~ — > f3 k . Since the partition function 
now does not factorize, as is the case for the GGE, the 
expectation values of products of conserved charges are 
not necessarily uncorrelated. From the partition function 
one can derive the probability Pjr (n) of having at least n 
particles in the state k. This is done by starting the sum- 
mation of rife at n in ( 13 1 and then using the recursion 
relation ( 14 ) to obtain 



Pc{n) = e- n P*Z N _ n /Z N . 



(15) 



The probability of having exactly n particles in the state 
k is therefore 



1 



P*(n) = — [e- n ^Z N - n -, 



-W>>Z N -^i) . (16) 



Similarly, one can derive the probability of having at least 
rik states in k and n q in q with k 7^ q: 

p k, q ( n k,n q ) = e- n ^»- n ^Z N _ {nk+nq) /Z N , (17) 

from which an expression for Pk,q{nk,n q ) can be ob- 
tained. Expectation values are now computed as follows: 



N 



(*I* fe )GCE = X)np fc (n). 



(18) 



n=l 



One can numerically solve these coupled equations in or- 
der to obtain values for /3fc. Once all fa are determined 
we can compute the other expectation values, for instance 



N N-m 

GCE = y2 V" n l n 2 pk,q{ni,n 2 ). (19) 
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FIG. 8. The long time average (LTA) of the g2{x) com- 
pared with the predictions of the generalized canonical en- 
semble (GCE), the generalized Gibbs ensemble (GGE) and 
the canonical ensemble (CE), for a system with N, L — 20 
and c = 1000. 



B. Results 

In fig. |8]we compare the LTA, GGE, GCE and CE for 
g%{x). As expected, the predictions from CE and GGE 
are completely off. The predictions of GCE agree ex- 
tremely well with the LTA see (fig. 7]), the relative error 
being less than 0.5% for all x. This may be ascribed 
to the fact that the sum-rule for g^ (x) is only saturated 
to 99.6%. The reason why the GCE works is not com- 
pletely obvious. We would like to stress that both for 
the GGE and the GCE only the expectation values of 
^l^fc are nxecl y i & corresponding Lagrange multipliers 
/3fc. In principle one could expand the GGE or GCE by 
including products of conserved charges k^ki^ k 2 

with their own Lagrange multiplier /?fc 1; fc 2 , but that is 
not what is done here. Apparently, by demanding that 
(A^ 2 ) = (A^) 2 one almost completely fixes higher order ex- 
pectation values. To gain more insight, we compare the 
expectation values k ^l +q ^ k+q ) for the LTA, GCE 
and the GGE plotted in fig. |9] for k = 0,2n/L. If we 
first focus on the top graph (k — 0), we see that for all 
q the LTA and GCE agree extremely well. The GGE is 
only valid for large q > 10, from which we conclude that 
at this point the conserved charges become uncorrelated. 
Furthermore, from the inset we see that the fluctuations 
of ^q^o differ significantly for the GGE. In the bottom 
graph (k = 2ir/L), the correlations are much weaker com- 
pared to the case k = 0. We see that the LTA and GCE 
differ slightly now, but the GCE is still better than the 
GGE. 

The results here are presented for finite size. In the 
thermodynamic limit one might expect that the predic- 
tion of the GCE and GGE coincide for local observables. 
On the other hand, correlations like (^Q^o^^g) for 
small q cannot be approximated using Wick's theorem; 
for these correlations the GGE remains invalid. To study 
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FIG. 9. Comparison of (^l^ k^t^^ h+q) for k = (top) and 
k — 2n/L (bottom) as a function of q for the LTA, GCE and 
GGE. N, L = 20 and c = 1000. 




FIG. 10. The predictions of the GCE, GGE and CE for 
N.L = 150 and c = 1000. 



how correlation functions like g2 (x) behave in the ther- 
modynamic limit, we compared the predictions of GCE 
and GGE for a system with N,L — 150 in fig. [lOj We 
see that the difference between GCE and GGE is still 
present for small x, but it is considerably smaller than 
for the case of N, L — 20. This leads us to believe that in 
the thermodynamic limit the GCE and GGE yield equiv- 
alent predictions for g<i (x). 
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IV. CONCLUSIONS AND DISCUSSION 

In this paper we have considered an interaction quench 
by turning off the interactions in the Lieb-Liniger model. 
By using the Heisenberg picture explicitly, the time evo- 
lution was computed via the correlation functions on the 
fully interacting initial state. The exact method we have 
used allowed us to study various correlation functions af- 
ter the quench at all time and length scales. Although 
this paper focused on the results for the Lieb-Liniger 
model, we would like to stress that the methods presented 
in this paper can be applied to any model for which the 
Heisenberg picture can be used to efficiently obtain the 
time evolution and for which the correlation functions 
on the initial state can be computed exactly. For exam- 
ple, one could prepare the system in the ground state of 
the anisotropic Heisenberg spin chain. By switching off 
the anisotropy term, the time evolution is accessible from 
that of free fermions. 

In case of the Lieb-Liniger model, we have studied sev- 
eral types of correlation functions after the quench. As 
initial state, the ground state of the Lieb-Liniger model 
was used for various interaction strengths c > 0. The 
time behavior of the correlation functions was studied 
both in the short and long time regimes. As expected, 
the results are most pronounced for large c, although the 
qualitative features do not strongly depend on the initial 
interaction strength. As a byproduct we have computed 
the four-point function of the Lieb-Liniger model for the 
first time. 

The long time average has been compared with vari- 
ous statistical ensembles. As expected, the canonical en- 
semble fails to make correct predictions, as is usually the 
case for integrable models. The GGE gives better results. 
However, for an intermediate number of particles it still 
fails to capture all the features of the LTA. This discrep- 
ancy can be understood by the fact that the fluctuation of 
the total number of particles is large, while for the quench 
under consideration the total number of particles remains 
constant. By introducing the GCE, which keeps the total 
number of particles explicitly fixed, correct predictions 
for the LTA are obtained. The failure of the GGE can 
also be explained by the fact that correlations between 
the conserved charges are not kept. The GCE, which 
is essentially the GGE with one additional constraint, 
seems to give the correct correlations between the con- 
served charges; why this is the case remains unclear on a 
more formal level. It is worth to mention that the GCE 
can be applied in the same cases as the GGE. A firmer 
test of the GCE would be to study even higher order cor- 
relation functions such as g n (x) = (("J/' (x)) n ('5 (0))") . In 



the thermodynamic limit it is expected that the particle 
number fluctuations of the GGE become irrelevant for lo- 
cal observables such as the 172 (z); a comparison of GGE 
and GCE for a large number of particles confirm this. 

It would also be interesting to see what the effect is of 
considering a different basis of conserved charges. For ex- 
ample in [55] the validity of the GGE was studied using 
two different bases for the conserved charges, of which 
only one agreed with the long time average. To make a 
connection with the interacting Lieb-Liniger model with 
c > 0, an appropriate basis would consist of the con- 
served charges in terms of derivatives of the transfer ma- 
trix, while sending c — > 0. This leads to the following set 
of conserved charges: Q n = J^ fe k n W k ^ k . However, for- 
tius quench problem the expectation values (4'\Qn\4 > } f° r 
n > 3 are not well-defined because of the 1/fc 4 behavior 
of (</>|\f r t^/.|<?!>) for large k [37] which reflects the ultra- 
locality of the Lieb-Liniger model. This does not imply 
that the GGE fails in this case, as one could regularize 
the results by introducing a lattice spacing a while keep- 
ing the system integrable. One can then write down the 
expectation values as function of a and send a — > at 
the very end of the calculation. This will be investigated 
in future work. 

A way to generalize the results of this paper is by 
considering quenches starting from arbitrary interaction 
strength c\ and ending in a different arbitrary ci. In 
this case, one cannot use the Heisenberg picture to com- 
pute the time evolution of observables as was done in 
this paper. One could try to solve the time evolution in 
the Schrodinger picture by making a spectral decomposi- 
tion of the initial state in terms of the eigenstates of the 
Hamiltonian after the quench. The first problem in this 
approach is the need for the overlap coefficients between 
the initial state and final states. For the Lieb-Liniger 
model, these overlap coefficients are only known in one 
particular case [SJ; the general problem is still unsolved. 
Secondly, the spectral decomposition typically involves 
an exponential number of states as function of the sys- 
tem size. Unless there are huge degeneracies present in 
the spectrum, as is only the case for the special points 
c = and c = 00, performing the spectral sum seems 
intractable. We will return to these and further issues in 
future publications. 

V. ACKNOWLEDGEMENTS 

This work is part of the research programme of the 
Foundation for Fundamental Research on Matter (FOM) , 
which is part of the Netherlands Organisation for Scien- 
tific Research (NWO). 



[1] M. Greiner, O. Mandel, T. Hansen, and I. Bloch, Nature [2] T. Kinoshita, T. Wenger, and D. S. Weiss, Nature 440, 
419, 51 (2002). 900 (2006). 



9 



S. Hofferberth, I. Lesanovsky, B. Fischer, T. Schumm, 

and J. Schmiedmayer, Nature 449, 324 (2007). 

A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalat- 

tore, Rev. Mod. Phys. 83, 863 (2011). 

T. Kinoshita, T. Wenger, and D. S. Weiss, Science 305, 

1125 (2004). 

A. H. van Amerongen, J. J. P. van Es, P. Wicke, K. V. 
Kheruntsyan, and N. J. van Druten, Phys. Rev. Lett. 
100, 090402 (2008). 

A. Imambekov, I. E. Mazets, D. S. Petrov, V. Gritsev, 
S. Manz, S. Hofferberth, T. Schumm, E. Demler, and 
J. Schmiedmayer, Phys. Rev. A 80, 033604 (2009). 
S. Manz, R. Bxicker, T. Betz, C. Roller, S. Hofferberth, 
I. E. Mazets, A. Imambekov, E. Demler, A. Perrin, 
J. Schmiedmayer, et al., Phys. Rev. A 81, 031610 (2010). 
V. Gritsev, T. Rostunov, and E. Demler, J. Stat. Mech. 
2010, P05012 (2010). 

D. Muth, B. Schmidt, and M. Fleischhauer, New J. Phys. 
12, 083065 (2010). 

A. Lamacraft, Phys. Rev. A 84, 043632 (2011). 

S. Sotiriadis, D. Fioretto, and G. Mussardo, J. Stat. 

Mech. 2012, P02017 (2012). 

E. H. Lieb and W. Liniger, Phys. Rev. 130, 1605 (1963). 
M. Girardeau, J. Math. Phys. 1, 516 (1960). 

P. Calabrese and J. Cardy, Phys. Rev. Lett. 96, 136801 

(2006) . 

P. Calabrese and J. Cardy, J. Stat. Mech. 2007, P06008 

(2007) . 

J.-S. Caux and P. Calabrese, Phys. Rev. A 74, 031605 
(2006). 

P. Deuar, A. G. Sykes, D. M. Gangardt, M. J. Davis, 
P. D. Drummond, and K. V. Kheruntsyan, Phys. Rev. A 
79, 043619 (2009). 

V. E. Korepin, N. M. Bogoliubov, and A. G. Izer- 
gin, Quantum Inverse Scattering Method and Correlation 
Functions (Cambridge, 1993). 

J.-S. Caux, P. Calabrese, and N. A. Slavnov, J. Stat. 
Mech. 2007, P01008 (2007). 

T. Kojima, V. E. Korepin, and N. A. Slavnov, Commun. 

Math. Phys. 188, 657 (1997). 

J.-S. Caux, J. Math. Phys. 50, 095214 (2009). 

M. Rigol, V. Dunjko, V. Yurovsky, and M. Olshanii, 

Phys. Rev. Lett. 98, 050405 (2007). 

E. T. Jaynes, Phys. Rev. 106, 620 (1957). 

M. A. Cazalilla, Phys. Rev. Lett. 97, 156403 (2006). 

M. Cramer, C. M. Dawson, J. Eisert, and T. J. Osborne, 

Phys. Rev. Lett. 100, 030602 (2008). 

T. Barthel and U. Schollwock, Phys. Rev. Lett. 100, 

100601 (2008). 

M. Kollar and M. Eckstein, Phys. Rev. A 78, 013626 

(2008) . 

D. Fioretto and G. Mussardo, New J. Phys. 12, 055015 
(2010). 

A. C. Cassidy, C. W. Clark, and M. Rigol, Phys. Rev. 
Lett. 106, 140405 (2011). 

P. Calabrese, F. H. L. Essler, and M. Fagotti, Phys. Rev. 
Lett. 106, 227203 (2011). 

D. M. Gangardt and M. Pustilnik, Phys. Rev. A 77, 
041604 (2008). 

T. Caneva, E. Canovi, D. Rossini, G. E. Santoro, and 
A. Silva, J. Stat. Mech. 2011, P07015 (2011). 
P. T. Landsberg, Thermodynamics - with quantum sta- 
tistical illustrations (Interscience Publishers, 1961). 



[35] P. Borrmann and G. Franke, J. Chem. Phys. 98, 2484 
(1993). 

[36] M. Wilkens and C. Weiss, J. Mod. Optic 44, 1801 (1997). 
[37] M. Olshanii and V. Dunjko, Phys. Rev. Lett. 91, 090401 
(2003). 



